Evolution of neutron stars with toroidal magnetic fields: Axisymmetric simulation in 

full general relativity 
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We study the stability of neutron stars with toroidal magnetic fields by magnetohydrodynamic 
simulation in full general relativity under assumption of axial symmetry. Nonrotating and rigidly 
rotating neutron stars are prepared for a variety of magnetic field configuration. For modeling the 
neutron stars, the polytropic equation of state with the adiabatic index F = 2 is used for simplicity. It 
is found that nonrotating neutron stars are dynamically unstable for the case that toroidal magnetic 
field strength varies cx with k > 2 (here ro is the cylindrical radius), whereas for k = 1 the 

neutron stars are stable. After the onset of the instability, unstable modes grow approximately in 
the Alfven time scale and, as a result, a convective motion is excited to change the magnetic field 
profile until a new state, which is stable against axisymmetric perturbation, is reached. We also find 
that rotation plays a role in stabilization, although the instability still sets in in the Alfven time 
scale when the ratio of magnetic energy to rotational kinetic energy is larger than a critical value 
~ 0.2. Implication for the evolution of magnetized protoneutron stars is discussed. 

PACS numbers: 04.25.Dm, 04.40.Nr, 47.75.+f, 95.30.Qd 
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I. INTRODUCTION 

Neutron stars observed in nature are magnetized with 
the typical magnetic field strength ~ 10^^-10"'^^ G 
The field strength is often much larger than the canoni- 
cal value as ~ 10^^ G for a special class of the neutron 
stars such as magnetars The field strength at the 
birth of neutron stars may be also much larger than the 
canonical value, because in the supernova gravitational 
collapse, rapid and differential rotation of the collapsing 
core could amplify the magnetic field. In the presence of 
a radial magnetic field , the toroidal field is am- 
plified by winding in the presence of differential rotation, 
and the field strength increases with time approximately 
according to (see, e.g., [1, 
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where we adopt the typical magnitude of the angular ve- 
locity and the typical cooling time of the protoneutron 
star for and t. The large field strength of the magne- 
tars may be generated by such process and subsequently 
be confined inside the neutron star for thousands of years 
0. This suggests that even for the normal pulsar, the 
toroidal field strength inside the neutron star may be 
much larger than the canonical value. Thus, strongly 
magnetized neutron stars may be common in nature. In 
particular, the toroidal field is likely to be much stronger 
than the poloidal fields inside neutron stars. In this pa- 
per, we focus on the effect of such strong toroidal mag- 
netic fields for the dynamical evolution of neutron stars. 

Stars with purely toroidal magnetic fields in a stably 
stratified structure are known to be unstable against the 
Tayler instability [SSBIl (see also ^pendix A). Ac- 
cording to a perturbative study in [1, 0, S, Q , the most 



unstable motions are driven by axisymmetric (to = 0) 
and nonaxisymmetric to = 1 modes with nearly horizon- 
tal displacement. The unstable modes are predicted to 
grow approximately on an Alfven time scale. The Alfven 
time scale of magnetized neutron stars are estimated to 
be very short as 
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where R and p are characteristic radius and density of 
the neutron star, va is the Alfven speed, and we use 
VA — -8-^/(477/9)1/2, rpj^g time scale for the growth of 
the Tayler instability is only by one order of magnitude 
longer than the dynamical time scale of neutron stars 
which is 1 ms. Thus, the instability associated with 
the strong magnetic fields may affect even early evolution 
of the protoneutron star and, consequently, supernova 
explosion. However, perturbative studies do not clarify 
anything in the nonlinear evolution stage reached after a 
sufficient growth of the instability. 

To understand roles of strong toroidal fields on the 
evolution of neutron stars and protoneutron stars, nu- 
merical simulation is probably the best approach. In this 
paper, we present our new numerical results obtained 
by general relativistic magnetohydrodynamic (GRMHD) 
simulation, for which our GRMHD code recently devel- 
oped [10] is used. We prepare neutron stars with purely 
toroidal magnetic fields in axisymmetric equilibria com- 
puted by a method described in ll]. As a first step to- 
ward a deep understanding of the Tayler instability, we 
focus on the to = mode imposing axial symmetry. As 
shown in Q, neutron stars are unstable if certain con- 
dition is satisfied for the magnetic field profile and for 
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the rotation rate. We confirm this fact in the present 
numerical simulation. In addition, we follow evolution of 
the unstable stars after the onset of the Tayler instability 
and show that associated with the growth of this instabil- 
ity, a convective motion is driven inside the neutron star. 
Then, the magnetic fields are redistributed and eventu- 
ally their profile relaxes to a new state which is stable 
against axisymmetric perturbation. 

The remainder of this paper is organized as follows. 
In Sec. II, we briefly review formulation and numeri- 
cal methods for our GRMHD simulations. Section III 
presents numerical results for nonrotating and rotating 
neutron stars separately. Section IV is devoted to a sum- 
mary and discussion about implication of the present re- 
sults on the evolution of neutron stars. In Appendix 
A, we present a result of linear perturbative study for 
neutron stars of purely toroidal magnetic flelds, which 
validates our numerical results qualitatively. Through- 
out this paper, we adopt geometrical units in which G = 
1 = c, where G and c denote the gravitational constant 
and speed of light, respectively. Cartesian coordinates 
are denoted by = {x,y,z). The coordinates are ori- 
ented so that the symmetric axis is along the ^-direction. 
We define the coordinate radius r = ^y x-^ + y"^ + z^, 
cylindrical radius zu = \/ x^ +2/^, and azimuthal angle 
93 = tan^^(?//x). Coordinate time is denoted by Greek 
indices /i, t^, • • • denote spacetime components, and small 
Latin indices • • • denote spatial components. 



[Tsj). the following dynamical gauge condition is employed 

dtOL = -aK^, (3) 
9t/3'=r^(F, +A<9tF,), (4) 

where a is the lapse function, /3' the shift vector, and At 
time step in numerical computation. 

A conservative shock-capturing scheme is employed to 
integrate the GRMHD equations. Specifically we use a 
high-resolution central scheme [l^ [l3| with the third- 
order piece-wise parabolic interpolation and with a steep 
min-mod limiter in which the limiter parameter 6 is set 
to be 2.5 (see appendix A of Multiple tests have 

been performed with these codes, including MHD shocks, 
MHD wave propagation, magnetized Bondi accretion, 
and magnetized accretion onto a neutron star (IQ] . This 
code has been already applied to the evolution of magne- 
tized hypermassive neutron stars to a black hole ^^.Tigj 
and to supernova gravitational collapse of strongly mag- 
netized and rotating core 4] , and derived reliable numer- 
ical results. 

In the present paper, we initially give a purely toroidal 
magnetic field. In such a case, poloidal magnetic fields 
are never generated in the axisymmetric spacetime. 
Thus, we only solve the toroidal field component. 

As initial conditions for the numerical simulation, we 
prepare magnetized neutron stars in equilibrium pH . For 
computing the equilibrium, we give the polytropic equa- 
tion of state as 



II. METHOD FOR NUMERICAL SIMULATION 
A. Formulation and methods 

The stability of magnetized neutron stars and the fate 
of unstable neutron stars are investigated by GRMHD 
simulation assuming that the ideal MHD condition holds. 
In this paper, we assume the axial symmetry and focus 
only on the Tayler instability against axisymmetric per- 
turbation. The simulation is performed by a GRMHD 
code for which the details are described in This 
code makes long-term numerical evolutions of relativistic 
magnetized neutron stars possible. It solves the Einstein- 
Maxwell-MHD system of coupled equations, both in axial 
symmetry and in 3-1-1 dimensions, without approxima- 
tion. The code evolves the spacetime metric using the 
BSSN formulation we evolve the conformal three 

metric 7ij = 7~^^^7ij, a conformal factor = ln(7)/12, a 
tracefree extrinsic curvature Aij — e^*'^(-ft'ij — 7^ Xj.'^/S), 
trace of the extrinsic curvature ^ and an auxiliary 
three variable Fi = '^jdj^fij- Here, jij is the three- 
metric and 7 = det(7ij). For axisymmetric simulation, 
the Cartoon method is employed [H, [13] : Namely, the 
Einstein equation is solved in the Cartesian coordinates 
imposing an axisymmetric boundary condition and the 
hydrodynamic equation is in the cylindrical coordinates. 

As in previous axisymmetric simulations (e.g., (lol . 



P - k/, (5) 

where P, p, k, and F are the pressure, rest-mass den- 
sity, polytropic constant, and adiabatic constant. In this 
work, we choose F = 2. Because k is arbitrarily chosen 
or else completely scaled out of the problem, we adopt 
the units of k = 1 in the following (i.e., the units of 
c = G = K = 1). 

In numerical simulation, we adopt the F-law equation 
of state 

P = (F - l)pe, (6) 
where e is the specific internal thermal energy. 

B. Diagnostics 

We monitor the total baryon rest mass M,, ADM 
(Arnowitt-Deser-Misner) mass M, and angular momen- 
tum J, which are computed, in axial symmetry, by 

M* = J pu^^/^(fx, (7) 
M = y pADuVld^^^ (8) 
J — phu*u^y^—g(fix, (9) 
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TABLE L List of characteristic quantities for neutron stars with toroidal magnetic fields. Value of k, central density, pc, baryon 
rest mass, M,, ADM mass, M, ratio of equatorial circumferential radius R to M, ratio of the rotational kinetic energy to 
the gravitational potential energy, Trot/W, ratio of the internal thermal energy to W, Eint/W, ratio of the electromagnetic 
energy to W, E-em/W , non-dimensional angular momentum parameter, J/M'^, central value of the lapse function, a^, angular 
velocity, Q, and Alfven time scale defined by Eq. (I21|l . All the quantities are shown in units of c = G = k = 1. For models 
AXY, BXY, and CXY, k — 1, 2, and 3, respectively. "X" denotes the value of lOpc and "Y"(=H[, M, L) denotes the relative 
strength of the magnetic field. Models RAXY and RBXY denote rapidly rotating neutron stars (meaning of A, B, X, and Y is 
the same as above) . Models MBXy denote moderately rapidly rotating neutron stars. Models EBXy and SBXy denote rotating 
models with very strong and strong magnetic fields, respectively, "y" approximately denotes lOOT/W . In the last column, the 
stability determined by the numerical simulation is described. 



Model 


k 


Pc 






M 


R/M 




Eint/W 


E-bm/W 










Q, 


ta/M Stable ? 


ASH 


1 


0.3000 


0.1793 


0.1637 


4.821 





0.5858 


0.0146 







0.4786 





67.8 


Yes 


A3L 


1 


0.3000 


0.1798 


0.1637 


4.749 





0.5931 


9.8 X 10" 


4 





0.4756 





256 


Yes 


A2H 


1 


0.2000 


0.1715 


0.1574 


5.599 





0.5148 


0.0135 







0.5732 





89.6 


Yes 


B3H 


2 


0.3000 


0.1794 


0.1637 


4.783 





0.5880 


0.0122 







0.4772 





73.4 


No 


B3M 


2 


0.3000 


0.1797 


0.1637 


4.750 





0.5927 


2.05 X 10" 


-3 





0.4757 





177 


No 


B3L 


2 


0.3000 


0.1798 


0.1637 


4.747 





0.5932 


9.3 X 10~ 


4 





0.4755 





263 


No 


B2H 


2 


0.2000 


0.1714 


0.1573 


5.543 





0.5168 


0.0108 







0.5715 





98.5 


No 


B2M 


2 


0.2000 


0.1717 


0.1573 


5.516 





0.5201 


3.11 X 10" 


.3 





0.5703 





182 


No 


B2L 


2 


0.2000 


0.1717 


0.1574 


5.505 





0.5210 


1.16 X 10" 







0.5700 





297 


No 


C3L 


3 


0.3000 


0.1798 


0.1638 


4.747 





0.5931 


1.16 X 10" 


-3 





0.4755 





236 


No 


C2H 


3 


0.2000 


0.1714 


0.1573 


5.537 





0.5169 


0.0108 







0.5713 





98.5 


No 


RA2H 


1 


0.2000 


0.1986 





1821 


6.732 


0.0808 


0.4555 


0.0145 




0.5667 





5383 


0.3159 


99.5 


Yes 


RA2L 


1 


0.2000 


0.2021 





1848 


6.491 


0.0866 


0.4577 


1.78 X 10" 


-3 


0.5894 





5318 


0.3271 


272 


Yes 


RB2S 


2 


0.2000 


0.1981 





1817 


6.283 


0.0796 


0.4559 


0.0176 




0.5614 





5382 


0.3144 


84.4 


No 


RB2H 


2 


0.2000 


0.2002 





1835 


6.594 


0.0839 


0.4549 


0.0126 




0.5783 





5352 


0.3211 


104 


Yes 


RB2L 


2 


0.2000 


0.2023 





1850 


6.478 


0.0869 


0.4575 


1.78 X 10" 


-3 


0.5906 





5315 


0.3276 


271 


Yes 


MB2H 


2 


0.2000 


0.1908 





1751 


5.863 


0.0611 


0.4708 


0.0163 




0.4870 





5460 


0.2852 


82.6 


No 


MB2M 


2 


0.2000 


0.1906 





1747 


5.797 


0.0597 


0.4744 


0.0108 




0.4817 





5456 


0.2832 


99.9 


No 


MB2L 


2 


0.2000 


0.1903 





1743 


5.738 


0.0581 


0.4780 


5.66 X 10" 


-3 


0.4758 





5453 


0.2809 


137 


Yes 


MB2L' 


2 


0.2000 


0.1856 





1700 


5.649 


0.0447 


0.4882 


4.66 X 10" 


-3 


0.4151 





5512 


0.2512 


149 


Yes 


EB27 


2 


0.200 


0.1915 





1765 


6.251 


0.0659 


0.4550 


0.0427 




0.5040 





5492 


0.2886 


54.8 


No 


EB25 


2 


0.200 


0.1858 





1713 


5.995 


0.0502 


0.4667 


0.0433 




0.4353 





5559 


0.2592 


52.6 


No 


EB23 


2 


0.2000 


0.1793 





1653 


5.798 


0.0302 


0.4833 


0.0392 




0.3343 





5635 


0.2078 


53.9 


No 


EB21 


2 


0.2000 


0.1745 





1610 


5.744 


0.0144 


0.4922 


0.0443 




0.2286 





5708 


0.1465 


50.6 


No 


SB24 


2 


0.2000 


0.1845 





1692 


5.702 


0.0430 


0.4854 


0.0137 




0.4052 





5536 


0.2459 


88.0 


No 


SB23 


2 


0.2000 


0.1805 





1656 


5.629 


0.0305 


0.4955 


0.0111 




0.3400 





5586 


0.2109 


97.2 


No 


SB21 


2 


0.2000 


0.1744 





1601 


5.572 


0.0108 


0.5091 


0.0116 




0.2005 





5671 


0.1286 


95.1 


No 



where is the four velocity, g is the determinant of 
the spacetime metric, h is the specific enthalpy {h = 

1 + e + P/p), and 



PADM = [ph{au'f ~ P]e-^ 



e ^ 
16^ 



Re 



-4<t> 



(10) 



Here, R is the Ricci scalar with respect to 7^ . Hereafter, 
initial ADM mass is denoted by Mq. 

In addition to the above quantities, we monitor the 
internal thermal energy -Eint, rotational kinetic energy 
Trot, total kinetic energy Tkin, and electromagnetic en- 
ergy T;em, written by 



Sir 



pu^e^ 



-gd X, 



Trot = 2 ' '^^^ U^pVt^/^d X, 



(11) 
(12) 



Tkin 



— / phu^Uiv''y/—gd'^x, 



-gd^x, 



(13) 
(14) 



where = h^h^^ is a magnetic vector in the frame co- 
moving with fluid elements (e.g., andw' = /u^. In 
this paper, magnetic field strength is defined by \/47r&2. 
We note that E-^m is defined originally by 



(15) 



where T^^ is the electromagnetic part of the energy mo- 
mentum tensor and is the hypersurface normal, and 
hence, the definition is different from that in [l^. This 
definition is based on mimicking the definition of -Eintj 
which is 



' hydro 
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where T^y^^^{= phu'^u" + Pgt''') is the non- 
electromagnetic part of the energy momentum tensor. 

Once each energy component is obtained, gravitational 
potential energy is defined by 

W = Ah + Eint + EeM + Tkin - M. (17) 

The Alfven speed in relativity is defined by 



62 



VA = 



ph + 62 ■ 



Associated Alfven time scale is 



TA 



L 

va' 



(18) 



(19) 



where i is a characteristic length scale and in the present 
context, it is approximately equal to stellar radius. Usu- 
ally, the Alfven time scale denotes a characteristic time 
during which Alfven waves propagate for the characteris- 
tic length scale. In the present case, we assume the axial 
symmetry and presence of purely toroidal magnetic field, 
and hence, Alfven waves play no role. Nevertheless, a 
dynamical instability analyzed in this paper grows in the 
time scale of order ta ■ For this reason, we here define an 
averaged Alfven time scale from global quantities as 



VA 



{ph + b'^)u*- X . 



int 



2E^ 



EM 



(20) 



where we use the relation h = 1 + Te which holds in the 
F-law equation of state. Then, we define averaged Alfven 
time scale as 

TA = ^, (21) 
VA 

where R is equatorial stellar radius. 



C. Initial condition 

We prepare a variety of neutron stars in equilib- 
rium changing the compactness, profile and strength of 
toroidal magnetic fields, and rotational kinetic energy. 
The initial conditions are derived in the same method as 
that described in [ll[. In the present case, we give the 
toroidal magnetic field according to the relation 



BQU^{pha^^^^)^ , 



(22) 



where k and Bq are constants which determine the field 
profile and field strength, respectively. Because of the 
regularity condition along the symmetric axis, k has to 
be a positive integer. 7^1^ is the (f(p component of 7^ and 



approximately proportional to tu^ near the symmetric 
axis. Because is a function of p, the magnetic field is 
confined inside the neutron star. 

In this work, we choose A; = 1, 2, and 3. Because 
hip is proportional to w'^^ near the symmetric axis, the 

toroidal field strength defined by B'^ = b^^^p^^"^ / is 
proportional to w'^^~^ . Namely, for small values of fc, the 
fields are confined near the symmetric axis. References 
S S IB (and also Appendix A) predict that stars with k = 
1 are stable against axisymmetric perturbation, whereas 
those with fc > 2 are unstable, although rotation could 
stabilize the unstable mode. 

Several key quantities which characterize the magne- 
tized neutron stars are listed in Table I. Bq is chosen 
so as to get 10^^ < E^m/W < 4 x 10"^ Yor typ- 
ical neutron stars of mass '-^ 1.4M0, VF ~ 6 x 10^"^ 
ergs. The electromagnetic energy is approximately writ- 
ten as i^EM ~ {B'^)'^R^ /i, and hence, the magnetic field 
strength we consider here is extremely large as 10^^-10^^ 
G for i? « 10 km. Such choice is done simply to save com- 
putational time (note that the time scale for the growth 
of unstable modes is proportional to {B'^)~^; see below). 
In all the cases, the magnetic field is strong but not strong 
enough to modify the stellar structure significantly; e.g., 
for the nonrotating case, the shape of the neutron stars 
is approximately spherical. 

Even from this extreme setting, we can derive a generic 
physical essence because scaling relation, associated with 
the magnetic field strength, holds for the evolution of 
the unstable neutron star. Namely, if the magnetic field 
strength becomes half, the growth time scale for the 
Tayler instability becomes approximately twice longer, 
although the qualitative properties about the evolution 
of the unstable star are essentially the same. Hence, the 
artificial choice of the large magnetic field is acceptable 
for deriving generic physical properties. 

Compactness of the neutron stars is determined from 
the conditions that the central density, pc, is 0.300 or 
0.200 (in units of k = 1). We note that the maximum rest 
mass and gravitational mass of spherical neutron stars for 
r = 2 are 0.1799 and 0.1637, and the corresponding cen- 
tral density is « 0.318. This implies that the nonrotating 
neutron stars with pc — 0.300 are close to the marginally 
stable point against gravitational collapse. Indeed, the 
rest mass and ADM mass for such models are close to 
the values of marginally stable stars (cf. Table I). We will 
show that our code can follow such extremely compact 
stars stably for a long time J> 3000Mo. By contrast, with 
Pc — 0.2, the ratio of the stellar radius to the ADM mass 
becomes ~ 5.5 for the nonrotating neutron stars, which 
is a typical magnitude for neutron stars (for a hypothet- 
ical value of ADM mass 1.35M0, i? « 11 km). Thus, 
in this paper, we consider very compact and reasonably 
compact neutron stars. 

We also prepare a variety of rotating neutron stars. In 
this paper, we focus only on rigidly rotating neutron stars 
with a moderate compactness; pc is fixed to be 0.200. For 
neutron stars with F = 2 polytrope, the maximum ratio 
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FIG. 1: Evolution of central density and central value of the lapse function for models A3H (solid curves), B3H (long-dashed 
curves), A2H (dashed curves), and B2H (dashed-dotted curves). The units of time is initial ADM mass Mo. 



of rotational kinetic energy to gravitational potential en- 
ergy, T/W J is ~ 0.09 for compact neutron stars with 
R/M K, 6 |20j. Here, at the maximum ratio, velocity at 
the equatorial surface of the star is equal to the Kepler 
velocity. Taking into account this fact, we prepare rotat- 
ing stars with < T/W ^ 0.09. Specifically, we consider 
the following four sequences for studying dependence of 
stability on the rotation rate. First we consider rapidly 
rotating stars with T/W = 0.08-0.09 or T/W w 0.06 
and with 0.001 <, Eem/W 0.02. The models in these 
categories are specified with models "R***" and "M***" 
(see the caption of Table I for the meaning of ***). In 
the third and fourth sequences, we approximately fix the 
values of Ebm/W as 0.04 and 0.01 but change the values 
of T/W for a wide range. The models in these categories 
are specified with model "E***" and "S***". By study- 
ing stability of these models, the dependence of stability 
criterion on T/W and Eem/W is clarified. 

III. NUMERICAL SIMULATION 

A. Choosing the grid points and atmosphere 

Numerical simulation was performed assuming the ax- 
ial symmetry as well as the equatorial {z = 0) plane sym- 
metry. For covering computational domain, a nonuni- 
form grid of the following grid structure is adopted for -cu 
and z: 

{iAx l<i< Nq 

iAx + £,AiAx (23) 
X log[cosh{(i - No)/Ai}] No<i<N. 

Here, x'^ denotes or z, Ax is the grid spacing in the 
inner region, and A^O: -^j a-nd ^ are constants which 
determine the grid structure of the outer part. In this 
grid setting, the inner domain with < zu < NqAx and 
< ^ < NqAx is covered by a uniform grid. Neu- 
tron stars are always covered in an inner region with 
r < 2NqAx/?>. Nq is chosen to be 150 in this paper. 



In the present simulation, mass ejection and expan- 
sion of stars do not occur in a remarkable manner, and 
hence, it is not necessary to resolve the outer region as 
accurately as the inner region where neutron stars are 
located. Thus, we prepare a rather large grid spacing for 
the outer region choosing Ai = 50 and ^ = 10. N is set to 
be 240. In the following, we present results with this grid 
setting for all the cases. With this setting, outer bound- 
aries along each axis are located at L « 800Ax which is 
approximately equal to eight stellar radii (« Si?) . These 
are large enough for excluding spurious effects from outer 
boundaries at least for t <^ 3000A/o- Indeed, we per- 
formed a simulation for N = 220 (i.e., L « 600 Ax) while 
fixing other parameters for the grid structure, and found 
that results depend very weakly on N (i.e., L; see Fig. 
ini). Nevertheless, for smaller values of L, the spurious ef- 
fects coming from the outer boundaries are serious: For 
L — SOOAa; « 3R and 600Aa; « 6R, the computation 
crashes eventually at t ^ 2000Mo and 2700Mo, respec- 
tively, in the chosen non uniform grid. However, note 
that the time at which computation crashes depends on 
the grid structure and for the uniform-grid case, the com- 
putation does not crash at t = 3500Mo even for L = 3R. 

For a convergence test, we performed additional simu- 
lations for model B3H, choosing uniform grid with N = 
240, 300, and 360. For each case, the equatorial radius 
of the neutron star is covered by 80, 100, and 120 grid 
points, respectively. Thus, for TV = 300, the grid resolu- 
tion for the inner region with i < 150 is the same as that 
for the nonuniform grid. In these uniform grids, outer 
boundaries are located at three stellar radii {L — 3R), 
but in these cases, the computation does not crash until 
t = 3500Afo. A simulation was also performed in the 
nonuniform grid with N — 220, in which L w 600Ax, as 
mentioned above. Comparison of the results for five sim- 
ulations indicates that the resolution of our typical grid 
setting and location of outer boundaries are fine enough 
to derive quantitative results within 20-30% error (see 
the last two paragraphs of Sec. IIIIB 1|) . 

Because any conservation scheme of hydrodynamics is 
unable to evolve a vacuum, we have to introduce an artifi- 
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FIG. 2: Snapshots for profiles of rest-mass density and mag- 
netic pressure, 6^/2, for model A3H which is dynamical stable. 



cial atmosphere outside neutron stars. We initially assign 
a small rest-mass density of magnitude pat = Pmax x 10^^ 
where pmax is the maximum rest-mass density of the neu- 
tron star. With such choice, the total amount of the rest 
mass of the atmosphere is about 10~^ of the rest mass 
of the neutron star. Thus, the accretion of the atmo- 
sphere onto the neutron star plays a negligible role for 
their evolution. 



B. Numerical results 

1. Nonrotating case 

We performed numerical simulations for all the models 
listed in Table I. All the simulations stably proceeded for 
a sufficiently long time to more than SOOOMq. For the 
case that unstable modes of MHD grow in the Alfven 
time scale, the simulation time is long enough to deter- 
mine the stability of each neutron star and to follow sub- 
sequent evolution after the onset of the instability for our 
chosen models. 

By the numerical simulations, we find that all the mod- 
els with k = 1 are stable and the profiles of density and 
magnetic field do not change significantly during evolu- 
tion. By contrast, all the nonrotating models with k > 2 
are unstable irrespective of the magnetic field strength. 
In this case, the magnetic fields are redistributed, and as 
a result, a convective motion is excited inside the neutron 
star. The unstable stars slightly expand and the central 
density decreases (cf. Fig. [T|). Eventually, the star re- 
laxes to a new state which is stable against axisymmetric 



perturbations. These conclusions hold irrespective of the 
magnetic field strength and compactness of the neutron 
stars. In the following, we describe characteristic features 
for stable and unstable stars showing numerical results 
for specific models. 

Figure [T] plots evolution of central density and central 
value of the lapse function for models ASH, B3II, A2II, 
and B2H. This illustrates that for models ASH and A2H, 
the neutron stars simply oscillate around their hypothet- 
ical equilibrium states. The oscillation is excited because 
the initial model deviates slightly from the true equilib- 
rium. The oscillation amplitude is larger for pc = 0.3 
than that for pc = 0.2. This reflects a fact that the neu- 
tron stars with pc — 0.3 are close to the marginally sta- 
ble point against gravitational collapse, and a small per- 
turbation induces a large deviation from the equilibrium 
state. In contrast to models ASH and A2H, for mod- 
els B3H and B2H, the central density and lapse quickly 
change at i ~ 500Mo, implying that the stars expand. 
This is due to the fact that the profile of magnetic fields 
is modified during the evolution. 

Figures [2] and [3] display snapshots for the profiles of 
density and magnetic pressure for models ASH and BSH 
at selected time slices. For model ASH in which the 
neutron star is stable, the profiles remain approximately 
static besides a slight oscillation. By contrast, model 
BSH is dynamically unstable against redistribution of 
magnetic fields: For t <^ 400Mo, the profile of the mag- 
netic pressure distribution gradually varies, and then, for 
400 ^ t/Mo ^ 1000, the magnetic fields are redistributed 
violently. As described in Appendix A, unstable modes 
grow near the equatorial plane as well as in a high lati- 
tude of z ^ Mo and ru ~ 2Mo. For t ^ lOOOAfo, the pro- 
file approaches to a new state which is stable against ax- 
isymmetric perturbation. The maximum magnetic pres- 
sure decreases after the onset of the instability. As a 
result, the pinching effect by the toroidal magnetic fields 
are weaken. This is the reason that the star slightly ex- 
pands and the central density decreases. 

Figure |4] plots square root of magnetic pressure, 
along cylindrical axis {w axis) at selected time 
slices {t/Mo = 0, 920, and 1849.1) for model BSH. For 
comparison, the profile for model ASH at t = is shown 
together. For model BSH in which k — 2, the magnetic 
field strength initially distributes approximately in pro- 
portional to zu^ for zu <, Mo- As a result of the growth 
of an instability, this profile changes, and eventually, the 
magnetic field strength becomes approximately propor- 
tional to w for w <, Mo- Indeed, in the relaxed state, the 
profile is similar to that for model ASH in which k — 1 
and \/Wj2 cx nj for uj < Mo- This indicates that mag- 
netized stars with k — 1 are attractors for the unstable 
star in axial symmetry. 

Figure [5] plots velocity vector fields {v^ ,v^) at selected 
time slices for model BSH. As shown in this figure, con- 
vective motion and circulation are excited during the 
nonlinear evolution after the onset of instability. The 
velocity of the convective motion becomes maximum dur- 
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FIG. 3: Snapshots for profiles of rest-mass density and magnetic pressure for model B3H. 



ing the nonlinear evolution of the instability at t ~ 600- 
700Afo and the maximum velocity of this motion is ^ 5% 
of the speed of light. It is interesting to point out that 
the convective motion is present even after the magnetic 
field profile approximately relaxes to a stable state. 

Figure [H] plots the evolution of the ADM mass, in- 
ternal thermal, kinetic, and electromagnetic energy for 
model B3H. This shows that the ADM mass is approx- 
imately constant, and internal thermal energy does not 
vary significantly. By contrast, the electromagnetic en- 
ergy decreases significantly after the onset of dynamical 
instability, and with the quick decrease, the kinetic en- 
ergy increases steeply. This is because the convective 
motion is induced by the dynamical instability. In other 
words, the electromagnetic energy is transformed into the 
kinetic energy. 

The kinetic energy increases up to 30% of the elec- 
tromagnetic energy for model B3II. (Possible error size 
of the kinetic energy is 20-30% as we discussed later; cf. 



Fig. 9.) This holds for models with pc = 0.3 and k — 2. 
For pc — 0.2 and k = 2, Tkin increases to ^ 0.5E-em 
and for models with k — S, Tkin increases to ~ O.GSem- 
All these results indicate that the kinetic energy of the 
convective motion can reach to a value approximately as 
large as the electromagnetic energy for the case that the 
instability grows. 

This result suggests that for a protoneutron star with 
strong toroidal magnetic fields (see Sec. I for discussion), 
the kinetic energy of the convective motion may reach 

Because the kinetic energy could reach to such a large 
value, the convective motion triggered by this type of in- 
stability inside a protoneutron star may affect supernova 
explosion. 

As mentioned above, this instability always sets in for 
k > 2 irrespective of magnetic field strength and com- 
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TABLE II: The growth time, r, of the Tayler instability for 
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FIG. 4: Profile of square root of magnetic pressure, 
along -cij-axis at selected time slices, t/AIo — 0, 920, and 1849.1 
for model B3H. For comparison, the profile for model A3H at 
t = is shown together. 



pactness of the neutron stars. However, the growth time 
scale of this instabiUty depends strongly on the magnetic 
field strength. Figure [7] plots electromagnetic energy as a 
function of time for several models with k = 1 (left) and 
k = 2 (right). For models A3H and A2H which are stable 
models, it remains approximately constant, whereas for 
models with fc = 2, it always decreases during the evolu- 
tion and the decrease time scale is shorter for the larger 
electromagnetic energy. 

To determine the growth rate of the dynamically un- 
stable mode, it is convenient to see Tkin because it is 
initially zero and increases purely due to excitation of 
the unstable modes. Figure [5] plots Tkin/M, as a func- 
tion of time for models B3H, B3M, and B3L (left) and 
for models B2H, B2M, and B2L (right). Note that for 
other unstable models, the behavior of the curve is qual- 
itatively the same. This figure shows that after the onset 
of the instability, the kinetic energy increases approxi- 
mately in an exponential manner with time as oc e*/'^ 
where r is a constant. This indicates that the instability 
is dynamical. This exponential growth holds for all the 
unstable models irrespective of the value of k and mag- 
netic field strength. Thus, we refer to r as the growth 
time scale. 

Table II lists approximate values of t for all the un- 
stable models. It is found that r increases systematically 
with the decrease of magnetic field strength. In Table 
I, we also describe the values of ta calculated using Eq. 
(PTjl . We find that the order of magnitude of r agrees 
with Ta well, and furthermore, the relation t/ta ~ 0.3- 
0.6 holds irrespective of the value of k and the magnetic 
field strength. This indicates that the dynamical insta- 
bility grows in the Alfven time scale. 

As we have reported, the instability sets in only for 
k > 2 and the growth time scale is approximately pro- 
portional to the Alfven time scale. These imply that this 



instability is indeed the Tayler instability [1, Here- 
after, we refer to this instability as the Tayler instability. 

In this paper, we input a very high magnetic field 
strength of ~ 10-^^-10^'' G. Because the scaling holds 
as shown above, the present result may be applied for 
neutron stars of canonical field strength ~ 10^^-10^^ G 
or magnetar field strength ^ 10^^-10^^ G. As shown in 
Eq. the Alfven time scale is ~ 10-100 ms for the 
magnetar field strength and ~ 1-10 s for the canonical 
field strength. Thus, the growth time scale of the Tayler 
instability is <^ 10 s for the field strength larger than 
~ 10^^ G. If this instability sets in for a neutron star, the 
electromagentic energy will be redistributed and trans- 
formed into kinetic energy in a short time scale (as short 
as or shorter than the rotation period). 

For illustrating that the results presented so far de- 
pend only weakly on grid resolution, we show numerical 
results with different grid resolutions and grid structures 
for model B3H. The chosen grid strcuture is described 
in Sec. HI A. Figure [H plots the evolution of the central 
density, the central value of the lapse function, the elec- 
tromagnetic, and kinetic energy, respectively. We also 
show violation of the Hamiltonian constraint and con- 
servation of the ADM mass. Here the definition of the 
violation of the Hamiltonian constraint is the same as 
that shown in Eq. (43) of [3l ; an averaged Hamiltonian 
constraint is defined by using the rest-mass density as a 
weight. We note that the ADM mass does not have to 
be conserved in axial symmetry, but in the present con- 
text with negligible gravitational radiation, it should be 
approximately conserved. 

We find that the numerical results are qualitatively 
the same irrespective of the grid setting, and also de- 
pend quantitatively weakly on it: At approximately the 
same time, the neutron star expands due to the nonlin- 
ear growth of the Tayler instability, and then, the density, 
the lapse function, and the electromagnetic energy relax 
to a stable state (besides oscillation around new equilib- 
rium values). The minumum electromagnetic energy and 
the maximum kinetic energy achieved depends weakly on 
the grid resolution and differences among five runs are 
at most 20%. For t ^ lOOOAfo, the violation of the 
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FIG. 8: Tkin/M. as a function of time for models B3H, B3M, and B3L (left) and for models B2H, B2M, and B2L (right). 



Hamiltonian constraint is larger for the case that non- 
uniform grid is employed, but this is purely due to the 
grid structure chosen [3J| and magnitude of the violation 
eventually relaxes to be small. These facts demonstrate 
that our choice for the grid resolution and the grid struc- 
ture is appropriate for studying this type of instability. 
One caution is that the kinetic energy in the late time de- 
pends strongly on the grid resolution. The likely reason 
is that with poorer grid resolutions, numerical viscosity 
dissipates the circulation, resulting in the suppress of the 
convective energy. This is illustrated by the fact that for 
the case of "360H", the kinetic energy is largest among 
five runs. Thus, in reality, the convective motion may 
be approximately constant for a time much longer than 
the Alfven time scale. However, accurately following the 
convective motion for a long time is not main subject of 
this paper, and hence, we do not touch this problem in 
detail. 



2. Rotating case 

Numerical simulations for rigidly rotating neutron 
stars were performed for all the models listed in Ta- 



ble I. We find that the stability criteria for the rotat- 
ing stars are different from those for nonrotating stars. 
As in the nonrotating case, stars with k = I are always 
stable against axisymmetric perturbation irrespective of 
magnetic field strength. Also, for many of stars with 
k = 2, the Tayler instability occurs and grows approxi- 
mately in the Alfven time scale. In contrast to the non- 
rotating case, however, stars with k > 2 may be stable 
for the rapidly rotating case; at least, for the first ~ 10 
Alfven time scale, we do not find evidence for occurring 
the Tayler instability. Figure [TO] plots evolution of elec- 
tromagnetic and convective kinetic energy (defined by 
Tkin - Trot) for models RA2H, RB2H, RB2S, EB21, and 
EB27. This shows that for the rapidly rotating models, 
RA2H and RB2H, the electromagnetic energy remains 
approximately constant and indicates that they are sta- 
ble, at least, in the time scale oi t — 2000A/o, more than 
Wta [l^. It is worth noting that for model RB2H, k — 2 
and electromagnetic energy is strong as Eem/W ~ 0.01. 
Nevertheless, it is a stable model: This illustrates that 
rapid rotation suppresses the Tayler instability. By con- 
trast, models RB2S and EB27, which have even larger 
electromagnetic energy Eem/W ~ 0.02-0.04, are unsta- 
ble even for a large value of Tro^/W ~ 0.07. 
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FIG. 9: Evolution of central density, central value of the lapse function, electromagnetic energy, kinetic energy, violation of 
the Hamiltonian constraint, and violation of conservation of the ADM mass for model B3H with five different grid structures. 
"2401", "2201", "240H", "300H", and "360H" imply nonuniform grid with iV = 240 and iV = 220, uniform grid with N = 240, 
300, and 360, respectively. The units of time is initial ADM mass A/q. Numerical results for Oc and pc are not shown for run 
2201 because they are approximately the same as those for run 2401 for t ^ 2500Mo. Computation for run 2201 crashes at 
t ~ 2700A'/o. 




FIG. 10: Sem/A/* and (Tkin -Trot)/A/* as functions of time for models RA2H, RB2H, RB2S, EB21, and EB27. For comparison, 
results for model B2H are plotted together. For models RA2H and RB2H, (Tkin — T'rot)/A#, remains to be smaller than 5 x 10~^. 



These results indicate that the stabihty is determined 
by the ratio i?EM/Trot; i-e., the Tayler instabihty grows in 
the Alfven time scale only when EEm/Trot is larger than a 
critical value ~ 0.2. To clarify this fact, we generate Fig. 
II II which summarizes the stability properties for rotating 
models with k = 2. This figure indeed suggests that the 
Tayler instabihty grows only for Eem/Tioi ^0.2 (note 
that the dashed line denotes EEu/Tiot — 0.16). This 
result is in qualitative agreement with the Newtonian 
analysis of Appendix A. 

As seen in Fig. [TUl the fraction of the decrease in 
the electromagnetic energy during evolution is smaller 



for larger value of T^-ct/W. Also, shown is that for larger 
ratio of T,-ot/ E-^m, the maximum value of convective ki- 
netic energy, Tkin — Jlot, is smaller and the convective 
motion damps more quickly for the rotating models. In 
particular, for the rapidly rotating case, the convective 
kinetic energy is smaller than the electromagnetic energy 
by two or three orders of magnitude. As mentioned in 
the previous section, the damping is partly due to the nu- 
merical dissipation. However, significant difference in the 
damping rates between nonrotating and rotating models 
suggests that the damping is primarily due to a rotational 
effect. All these facts indicate that rotation plays a role 
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FIG. 11: Results of stability properties for rotating models 
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and the instability grows in the time scale of ^ ta- For the 
stars below the dashed line, instability does not grow in the 
Alfven time scale, although we cannot exclude the possibility 
that the stars are unstable and the unstable modes grow for 
a time scale S> ta- 



for stabilizing the axisymmetric unstable mode. 

As discussed in [1, Q, the Tayler instabihty induces 
a motion perpendicular to the rotation axis, and stabi- 
lization by rapid rotation is due to the presence of the 
Coriohs force. The Coriohs force pushes a fluid element, 
which deviates from its equilibrium location to the pos- 
itive cylindrical-radial direction, to the counter-rotation 
direction. As a result, the centrifugal force of the fluid 
element is weaken, and it is enforced to return toward its 
equilibrium position. Thus, it is natural that the Coriolis 
force suppresses the onset of the Tayler instability as well 
as convective motion in the meridian plane. 



IV. SUMMARY 

In this paper, we have reported stability of neutron 
stars with purely toroidal magnetic flelds. The stability is 
determined by performing GRMHD simulation. For the 
simulation, we prepare a variety of equilibrium neutron 
stars changing their compactness, strength and proflle of 
toroidal magnetic fields, and rotational kinetic energy. 
The following is the summary of the numerical results. 

(i) Magnetized stars with k = 1 are stable against 
axisymmetric perturbation irrespective of the magnetic 
field strength and the rotational kinetic energy. 

(ii) For the nonrotating case with fc = 2, magnetized stars 
are dynamically unstable irrespective of the magnetic 
field strength, and the magnetic field is redistributed ap- 
proximately in the Alfven time scale. The resulting pro- 
file of the magnetic fields is similar to that for k = 1, 
indicating that stars with A; = 1 are attractors if only the 
axisymmetric perturbation is taken into account. 

(iii) During the growth of the dynamically unstable 



modes, electromagnetic energy is transported to the ki- 
netic energy via the Tayler instability, and as a result, a 
convective motion is excited. The magnitude of the con- 
vective kinetic energy becomes approximately as large as 
the electromagnetic energy at the time when the non- 
linear growth saturates, for not-rapidly rotating neutron 
stars. 

(iv) For the case that a neutron star is rapidly rotat- 
ing, the Tayler instability is suppressed and the star with 
k = 2 may be dynamically stable against axisymmetric 
perturbation at least in about ten Alfven time scale. The 
numerical results suggest that if the rotational kinetic 
energy is more than ~ 6 times larger than the electro- 
magnetic energy, the star is stabilized. 

(v) Even for the unstable rotating models, convective mo- 
tion is not induced as remarkablely as in the nonrotating 
case. For the rapidly rotating case in which the rotational 
kinetic energy is larger than the electromagnetic energy 
by a factor of ^ 2, the maximum convective energy is 
smaller than the electromagnetic energy by two or three 
orders of magnitude. This indicates that rotation in gen- 
eral plays a role in stabilizing the axisymmetric Tayler 
instability. 

As mentioned in Sec. 1, protoneutron stars are likely 
to have strong toroidal magnetic fields if the progeni- 
tor of supernova gravitational collapse is rotating. In- 
deed, a number of recent MHD simulations of super- 
nova collapse of magnetized rotating stars have shown 
that the toroidal magnetic field is amplified by many 
order of magnitude during collapse and during subse- 
uent relaxation stage of the formed protoneutron star 
^ Hi m, H, H, m, HE m m, m, . The key mech- 
anism in this amplification is transport of the rotational 
kinetic energy to the electromagnetic energy via wind- 
ing induced by differential rotation. As indicated in 3, 
the electromagnetic energy could be comparable to the 
rotational kinetic energy at the end of the amplification. 
Our present numerical experiment suggests that when 
the condition i^EM/Trot <^ 0.2 is achieved, the protoneu- 
tron star may be subject to the Tayler instability. 

Many of the MHD studies for supernova collapse focus 
on the amplification of magnetic pressure associated with 
the amplification of the toroidal magnetic field strength, 
which increases the pressure behind shock waves and can 
help supernova explosion or drive a strong outflow along 
the rotational axis. Our present study suggests that the 
strong toroidal magnetic field may play a role not only 
in increasing the pressure for pushing shock waves but 
also in exciting a convective motion through the Tayler 
instability. After the onset of this instability, a large frac- 
tion of electromagnetic energy may be transported into 
the convective kinetic energy. Then the convection may 
help to carry a hot material near the surface of a pro- 
toneutron star toward the gain region and push stalled 
shock waves outward [Slj. As Eq. p4)) indicates, the 
kinetic energy of the convective motion may increase to 
^ 10^° ergs if the toroidal field strength becomes 10^^ G 
and the electromagnetic energy becomes as large as the 
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rotational kinetic energy. This value of the kinetic en- 
ergy amounts to ~ 10% of the energy required to drive 
supernova explosion. As we show in Fig. [5l vorticity 
is generated associated with the convective motion. If 
this vorticity is dissipated by viscosity, a large thermal 
energy is also generated. Such thermal energy may also 
contribute to pushing stalled shock waves (see similar dis- 
cussion in [32]). All these possibilities suggest that the 
Tayler instability should be taken into account in magne- 
torotational explosion scenarios for supernova explosion. 

In reality, electromagnetic energy in protoneutron star 
at birth is likely to be much smaller than rotational 
kinetic energy. Thus, at its birth, the protoneutron 
star will be stable against the Tayler instability. Sub- 
sequently, the toroidal magnetic fields are amplified by 
winding caused by differential rotation, and as a result, 
the electromagnetic energy will reach a magnitude com- 
parable to the rotational kinetic energy. Then, the pro- 
toneutron star could be unstable against the Tayler insta- 
bility. This suggests that this instability may play a role 
for stopping the growth of the toroidal field strength. At 
the end of the toroidal field amplification, the resulting 
electromagnetic energy is likely to be at most comparable 
to the rotational kinetic energy, and hence, the instabil- 
ity may not be as strong as that in the nonrotating stars, 
as illustrated in Sec. Ill B. Also, the convective energy 
driven by the instability is likely to depend strongly on 
the amplification process of the toroidal magnetic fields 
via winding. For example, if the degree of differential 
rotation is largest near the rotation axis, the instability 
will not be strong because the configuration is likely to 
be similar to that oi k = 1. On the other hand, the 
degree of differential rotation is largest at an interior of 
a protoneutron star far from the rotation axis, the Tay- 
lor instability will occur. The efficiency of the winding 
depends strongly on the magnetic field profile and the ro- 
tational profile of progenitor. To understand the role of 
the Tayler instability, well-resolved, long-term MHD sim- 
ulations of stellar core collapse have to be systematically 
performed for a number of initial conditions. 

Neutron stars, observed as ordinary pulsars, typically 
have rotation period 0.1-1 s and magnetic field strength 
of 10^^-10^'^ G. This class of neutron stars are not sub- 
ject to the Tayler instability. This is because the ratio 
of electromagnetic energy to rotational kinetic energy is 
much smaller than unity as 
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where / is moment of inertia. 

By contrast, magnetars of rotation period 5-12 s and 
of magnetic field strength lOi'^-lO^^ G [1| are subject 



to the Tayler instability. Present numerical results im- 
ply that stable magnetars should have a magnetic field 
profile which is at least stable against axisymmetric per- 
turbation. 

In this paper, we focus only on the stability against 
axisymmetric perturbations. This work should be re- 
garded as the first step toward deeper understanding of 
the Tayler instability in magnetized neutron stars. As 
shown in [1, 0, 0, Q , neutron stars with toroidal magnetic 
fields are also unstable against nonaxisymmetric pertur- 
bation. Nonaxisymmetric Tayler instability may grow in 
a different manner from axisymmetric one, and hence, 
the dynamical evolution process of neutron stars during 
the growth of the unstable modes as well as the final fate 
could also be different. Furthermore, this instability will 
occur even for k = \. In the axisymmetric simulation, 
only the stars with k > 2 are unstable and magnetic 
field profile of such unstable star eventually relaxes to a 
profile similar to that with k = 1. However, in three di- 
mensions, such star will be still unstable. This suggests 
that stars may never reach a stationary state. To answer 
this question, GRMHD simulation in full three dimen- 
sions is required. We plan to perform three dimensional 
simulation in the next step. 
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APPENDIX A: PERTURBATIVE ANALYSIS 

In this section, we present a result of perturbative anal- 
ysis on criteria for the onset of axisymmetric instabilities 
of neutron stars with toroidal magnetic fields. For the 
analysis, we follow Ref. [tJ. Newtonian gravity is as- 
sumed for the sake of clarity and simplicity, and thus, 
our purpose is to derive an approximate criterion for the 
instabihties. 

Basic equations describing ideal MHD are given by 



(Al) 



1 

+ — if-- V J 
47r ■' 

dtB' ^Vj{v'B' ~v^B'), 
y^B' = , 



B^VjB, - pg* , (A2) 

(A3) 
(A4) 



where p denotes the rest-mass density, the fluid veloc- 
ity, P the pressure, B^ the magnetic field, g* the gravita- 
tional acceleration, and the co variant derivative with 
respect to cc'. 
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We derive linear perturbation equations for rigidly ro- 
tating stars with purely toroidal magnetic fields in equi- 
librium. The velocity and the magnetic field for the equi- 
librium stars are written as 



(A5) 
(A6) 



where f2 and B{tu, z) are the angular velocity and the 
magnetic field strength, respectively. Here, we used the 
cylindrical polar coordinates (137, if, z). 

In order for the stability analysis to be tractable, we 
only consider axisymmetric perturbations of very short 
wavelength both in the w and z directions. Here, the 
short wavelength implies that the wavelength. A, of an 
oscillation mode is smaller than 5vt where 5v is a typical 
magnitude of the perturbed velocity field and r is a typi- 
cal change time scale of stellar structure. In other words, 
we perform the local analysis. We also employ the Cowl- 
ing approximation, in which perturbations of the gravity 
are omitted. 

In the short-wavelength approximation, linear pertur- 
bation equation for the mass conservation equation (|Aip 
is 



(A7) 



where 6Q denotes the Euler perturbation of the physi- 
cal quantity Q and we assume that ^ \dv'^ /m\ 
because of the short-wavelength approximation. Equa- 
tion (|A7p implies that the effect of density perturbation 
does not play a role. This is also because of the short- 
wavelength approximation imposed. Thus, sound waves 
or p- modes are filtered out in this analysis. 

In the short-wavelength approximation, the vj and z 
components of Eq. (|A2p become the following same equa- 
tion 



5P+—WB SB'-^ = . 
47r 



(A8) 



The other pieces of independent information extracted 
from Eq. (|A2[) are given by 

dtSv,^ + 2vjn = -5- SB^dj {wB) , (A9) 
Ait 

p[dt{dz5v^ - d^dvz) - 2-070. d^Sv'^] 

= -L (_2S d.SB"^) - g^d.Sp + g.d^dp ,(A10) 
47r 

where gi is the apparent gravity, defined by 

g, EE g* + u^WjU, = (5; - ojn^, 0, <?:) . (All) 
The induction equation (|A3[) gives 



dtSB"^ ^ , dtSB' = , 



(A12) 



dtSB"^ = — dt6p - pSv^d. ( —] . (A13) 
rup X'^pJ 



In this study, we focus on adiabatic oscillations. Then, 
the relationship between SP and Sp is 



6P+e^^p = — {5p + e^^p), 

p 

where F denotes the adiabatic index, defined by 
'lnP\ 



r = 



(A14) 



(A15) 



and the Lagrangian displacement, which obeys 

5v' = dtC + v^djC - £.^djv' . (A16) 

For the axisymmetric perturbation with = QSJ, Eq. 
(|A16|) reduces to 



Then, Eq. (|A14p becomes 

-1- dtSP = - dtSp + Sv'A, 
vip P 

where Vs is the adiabatic sound speed, defined by 

'PF' 



(A17) 



(A18) 



(A19) 



and Ai the Schwarzschild discriminant. 



1 



Ai = c)i In p - - 9i In P . 



(A20) 



As shown in Eq. \kl2\ . we have 5B^ = SB"" 
because we are not interested in time-independent per- 
turbations. Thus, Eqs. (|X7l) - Cno| . (|XT3)) . and (IXTS)) 
are six independent equations for the six independent 
variables 5B'^ , 6p, and SP. In the local analysis, 
axisymmetric perturbations can be written as 

SQ{t, tu, z) — Qo exp{i{—at + Im + nz)} , (A21) 

where Qo is a constant, a the oscillation frequency, and 
(l,n) the meridional wavenumber vector. 

Substituting Eq. (|A2ip into the perturbation equa- 
tions, we obtain the following dispersion relation for a: 



1 + I -J ^ = 4^^-" 1 



^ ' Ah 



m J vj dq \wp; 



Here, va denotes the Alfven speed 
s the total meridional wavenumber 



(A23) 



(A24) 
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the apparent gravity along the constant phase or crests 



I 

9 = 9tx, ffz , 

n 



(A25) 



Ah the Schwarzschild discriminant along the constant 
phase or crests 



Ah = A^ , 



l_ 

n 



(A26) 



and d/dq the derivative along the constant phase or 
crests 



d _ d Id 
dq dm n dz 



(A27) 



The first, the second, and the third terms in the right- 
hand side of Eq. (|A22p are related to effects of the rigid 
rotation, the stratification (the buoyancy), and the mag- 
netic buoyancy, respectively. Here, it should be empha- 
sized that the effect of rigid rotation operates as a stabi- 
lizing agent because the first term in the right-hand side 
of Eq. (IA22P is always positive. 

The stabilities are determined by the dispersion rela- 
tion (|A22p . Specifically, the sign of its right-hand side 
determines the local stabilities; the stars are locally sta- 
ble (u nstable) if ct^ > q (0-2 < q). From Eqs . jM !])- 
([MT]) . we find that the right-hand side of Eq. (|M2l) is a 
quadratic in l/n, given by 



1 + 



( - I +c, 

n 



where 



v1 d 



z 



zup 



(A28) 



(A29) 



2vl 

Az + dzA^ 



^ In 



w? dz 



vjp 



v1 dw V'^P / ' 



(A30) 



4f]2 1 



^ ' Ar. 



^Inr^KASl) 



w? dz 



zup 



We therefore see that the stability condition <t^ > for 
any value of l/n is equivalent to the condition that the 
three inequalities a > 0, c > 0, and b^—Aac < are simul- 
taneously satisfied. Contrapositively, it is found that the 
star is unstable if any of the following three conditions is 
satisfied: 



4n'{l + ^]-{g^ + ^]A^ 



2^2 

ZD 



9z 



2«?\ vl d 



9z <: -A 



ru ) vl dzu 
^ B 
Vg dz \z17p 



—] <0,(A32) 
zup J 



d , 
In 



< 0. 



(A33) 



2 ) A 



dz 



zup 



dz \Zjop J 



Vg dz \zup 



A. 



< .(A34) 



Here, the final equation (|A34|) is equivalent to the condi- 
tion ^2 - 4ac > 0. 

Neutron stars are likely to be stably stratified because 
of their strong composition gradient [33[. As a result, the 
buoyancy inside the neutron star exerts as a stabilizing 
force as shown in Eqs. (|A32p and (|A33p . Equation (|A32p 
also shows that the criterion of the magnetic instability 
depends on whether the region considered is located in- 
side or outside the critical surface whose cylindrical ra- 
dius is defined by 



5ro 



(A35) 



Inside (Outside) the critical surface, if ^ In {B/zup) > 
{-^ln{B/zup) < 0), the third term of Eq. (|X32)) be- 
comes negative and the instability is promoted. 

Finally, we examine the magnetic stability of a par- 
ticular model, a slowly rotating star containing toroidal 
magnetic fields. We assume that matter distribution of 
the star is spherical (namely the magnetic and centrifu- 
gal forces are not strong enough to modify this spherical 
shape). As an example, an n = 1 polytropic sphere is 
considered because the analytic solution is available. Its 
matter distribution is given by 



P = Po 



P = Po 



sm{r/ro) 

sin(r/ro) 
(r/ro) 



(A36) 
(A37) 



where pa and Pq are the central values of the density and 
the pressure, respectively, and tq is stellar radius 



ro 



2Po 
A-kGpI 



For the magnetic field distribution, we take a simple 
form, as in [ll|, as 



B = 6o(p/po)'(r-sin0/ro) 



2/0-1 



(A38) 



where 60 and k are constants. Regularity of the mag- 
netic fields around the magnetic axis requires that k is 
a positive integer. Here, we used the polar coordinates 
(r, 0, (fi). Note that this magnetic field distribution is the 
Newtonian limit of that used in the present GRMHD sim- 
ulation. In this analysis, we omit the buoyancy or take 
= to focus on the magnetic instability. 
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First, we pay attention to the nonrotating case. Then, 
the instabihty conditions (|A32p and (|A33[) are written as 



Di{r,e) 



_ '0 



X In 



vi 



X In 



rsmt 
B 

p r sin 9 

cos 6— 
or 

B 

p r sin 9 



sin 9 



d cos 9 d 



dr 



<0, 



sin 6' d 
~d9 



<0, 



(A39) 



(A40) 



Note that by definition, Di and D2 are dimensionless 
quantities and are independent of the magnetic field 
strength bo. The left-hand side of Eq. (|A34[) vanishes 
in the present situation and this third instability condi- 
tion cannot give any useful information. For the case 
oi k = 1, it can be seen that Di = = D2 because 
di \n{B/p r sin 9) oc k — 1. Thus, the star is neutrally sta- 
ble for A: = 1. In Fig. [121 we give the distributions of 
Di and D2 on the meridional cross section for the case 
of fc = 2. In this figure, the darker regions have locally 
larger growth rates of the unstable mode, whereas the 
white regions correspond to neutrally stable ones (regions 
of Di ~ and D2 ^ 0). We then see that the magnetized 
stars with k = 2 are indeed locally unstable. In Fig. [T2l 
we confirm that the unstable regions determined by the 
criterion Di are separated by the the critical surface. 

The instability determined by Di occurs primarily near 
the equatorial plane and relatively high-density region. 
By contrast, the instability determined by D2 occurs near 
the surface and for the region of relatively weak magnetic 
field, and this indicates that this mode plays a minor role 
for redistribution of the magnetic field and for inducing 
convection. The instability found in our numerical simu- 
lation is likely to be associated with the mode determined 
by Di. Hence, in the following, we focus primarily on this 
mode. 

The modes associated with Di and D2 determine the 
instability in particular for l/n and l/n — > 00, re- 
spectively. Thus, we focus on the case for small values 
of l/n. As discussed in d, the Tayler instability is 
associated with a motion perpendicular to the magnetic 
axis. This type of motion corresponds to the limit of 
l/n — > in this analysis because Sv^ = —{l/n)Sv^ (see 
Eq. (jA7[) ). Thus, from this point of view, it is reasonable 
to pay attention for the mode associated with Di. 

Because the Tayler instability occurs for the magnetic 
field with /c = 2 as discussed above, henceforth, we only 
consider the case of fc = 2. For this model, the averaged 
Alfven speed va is given by 



VA = 



An / pd^x 



(315- 



2800 
« 1.25 va,o. 

where va is defined as 



2007r2 + 327r^; 



VA,0 



(A41) 



VA,0 



bo 



\/47rpo 



The growth time r for the Tayler instability defined in 
Sec. Ill is associated with increase in the kinetic energy. 
The growth time r for the most unstable mode is there- 
fore given by 



t/ta = (2crfA)"^ ~ 0.199 Min 



VA,0 

rod 



(A42) 



where Min((5) denotes the minimum value of Q{r,9). 
In the weak magnetic field approximation assumed, 
VA,ai'^o cr)~^ is given by 



VA,0 

ro 



(A43) 



where b — (rQU^Q)fe. In Fig. [T31 we show the growth 
time t/ta obtained by the local analysis as a function of 
l/n. As mentioned above, we focus on the case that l/n 
is small. Then, Fig. [13] shows that the unstable mode 
characterized hy l/n — is the most unstable one, whose 
growth time is given by 



t/ta w 0.089 . 



(A44) 



Thus, the minimum growth time is by a factor of 3-5 
shorter than that obtained by the GRMHD simulations 
(compare with Table II), but the order of magnitude 
agrees. For modes with a moderate value oll/n, as shown 
in Fig. I13[ the growth time t/ta increases. For a mode 
with l/n w 2, for example, t/ta ~ 0.2, which is 1/3- 
2/3 of the growth time shown in Table II. This result is 
reasonable because we here assume a Newtonian model 
whereas in the simulations, we adopt highly general rela- 
tivistic neutron stars for which the profiles of density and 
magnetic field are significantly different from the Newto- 
nian ones. 

Next, we consider the slowly rotating model. In the 
slow rotation and weak magnetic field approximation, the 
criterion of the Tayler instability becomes 




Mi 

. n 



dA < 0(A45) 



For the fc = 2 model, the rotational kinetic energy Trot 
and the electromagnetic energy E-em are, respectively, 
given by 



T - ^ 

J rot — :r 



pr 



2 sin^ 9n^ d^x 
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FIG. 12: Contours curves (solid curves) of D\ (left) and D2 (right) on the meridional cross section for the n = 1 polytropic star 
containing weak toroidal magnetic fields with k = 2. The darker regions are locally more unstable, whereas the white regions 
correspond to neutrally stable ones (regions with Di « « 02)- The contours of equi-D are linearly spaced; the difference 
between two adjacent contours is a sixth of the difference between the maximum and minimum values of D\ and D2- The 
maximum value is zero and the minimum values are —4.96 and —1.07 for Di and D2, respectively. The thick quarter circle 
denotes the surface of the star. Inside the star, the dashed curves and the dashed-dotted curve show equi-B contours and the 
critical surface, respectively. 
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FIG. 13: Growth time t/ta for the most unstable mode, given 
as a function of l/n. 




10 15 



FIG. 14: Critical ratio of the electromagnetic energy to the 
rotational kinetic energy, {Eem/Tioi)^, given as a function of 
l/n. 



-TT (TT 



E 



EM 



50.9 Poring 



(A46) 



1 

8^ 

5600 
2.45 rlbl 



(315 - 2007r^ + 327r*)r^6S 



(A47) 



The ratio of the electromagnetic energy to the rotational 
kinetic energy is then written as 



^EM/Trot ~ 0.605 



(A48) 



In terms of £'EM/7rot, thus, Eq. (|A45|) is rewritten as 

-EEM/Trot > (i?EM/Trot)c 



Min 



-2.42 



(A49) 



For the polytropic models with n — 1 and k — 2, numer- 
ical values of (-EEM/7rot)c ^'''^ shown as a function of l/n 
in Fig. [H 

Form this figure, it is found that for the most unstable 
mode, whose value of l/n is zero, the Tayler instability 
sets in if the condition. 



EEu/Trot > 0.49 , 



(A50) 



is satisfied. For modes with a moderate value of l/n, 
values of {EEM/Trot)c decreases, e.g. (£^EM/7'rot)c ~ 0.21 
for a mode with l/n 2. As argued in Sec. Ill, the 
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instability condition obtained by the GRMHD simulation 
is i?EM/7iot ^ 0.2, which is the same order as that of Eq. 
(IA501) . 

For the small values of Eem/Tioi, the unstable modes 
should have larger values oil/n, and correspondingly, the 
growth time scale becomes longer. This suggests that 
even if a model star appears to be stable in a numerical 
simulation for a finite duration, the star might become 
unstable for a sufficiently long run. 



As found from Fig. [Ml magnetized star is always un- 
stable for modes with \l/n\ « oo irrespective of the rota- 
tion rate. As mentioned previously, however, these modes 
are associated with I?2- This instability occurs near the 
stellar surface and its effect would not be significant for 
global redistribution of the magnetic field profile. Thus, 
although all the magnetized stars with k = 2 are unstable 
strictly speaking, this type of instability does not seem 
to play a significant role. 
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Here, let us explain why Hamiltonian constraint violation 
behaves like in Fig. 9(e). At t = 0, the degree of violation 
is much larger than 0.3 %; it is about 1%. The reason for 
this is that the initial condition does not satisfy the con- 
straint with a high accuracy because the initial condition 
is prepared by linear-interpolating the quantities of the 
equilibrium solution computed in a different grid struc- 
ture from that in the simulation. For t > 0, the degree of 
constraint violation decreases by the effect of constraint 
damping term included in this code. This term makes 
the equation for the conformal factor (</)), which primar- 
ily determines the Hamiltonian constraint violation, be 
parabolic. Therefore, the equation is not local one and 
the decrease rate of the constraint violation should de- 
pend on the grid structure. 

We cannot exclude a possibility that the rapidly rotating 
stars are eventually unstable after a longterm evolution. 
As shown in 0, [g|, the growth rate is proportional to 
v\/FI?Q, for the rapidly rotating case, and hence, it might 
be necessary to perform simulation for a much longer 
time scale. 

In this Appendix, we recover the gravitational constant 
G. Also, c does not denote the speed of light. 



